Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 2 February 2008 (MN MTpX style file v2.2) 



The Evolution of Binary Fractions in Globular Clusters 



N. Ivanova x * , K. Belczynski 2 f, J. M. Fregeau 1 , & F. A. Rasio 1 

1 Northwestern University, Dept of Physics & Astronomy, Evanston, IL 60208, USA 

2 Present address: New Mexico State University, Department of Astronomy, 1320 Frenger Mall, Las Cruces, New Mexico 88003-8001, USA 



in 
o 
o 

(N 
G 



> 

m 

7— I 

o 
in 
o 

Oh 

6 



X 



2 February 2008 



ABSTRACT 

We study the evolution of binary stars in globular clusters using a new Monte Carlo approach 
combining a population synthesis code (Star Track), and a simple treatment of dynami- 
cal interactions in the dense cluster core using a new tool for computing 3-body and 4-body 
interactions (Fewbody). We find that the combination of stellar evolution and dynamical 
interactions (binary-single and binary-binary) leads to a rapid depletion of the binary popu- 
lation in the cluster core. The maximum binary fraction today in the core of a typical dense 
cluster like 47 Tuc, assuming an initial binary fraction of 100%, is only about 5-10%. We 
show that this is in good agreement with recent HST observations of close binaries in the 
core of 47 Tuc, provided that a realistic distribution of binary periods is used to interpret the 
results. Our findings also have important consequences for the dynamical modeling of glob- 
ular clusters, suggesting that "realistic models" should incorporate much larger initial binary 
fractions than has usually been done in the past. 

Key words: binaries: close - binaries: general - methods: TV-body simulations - globular 
clusters: general - globular cluster: individual (NGC 104, 47 Tucanae) - stellar dynamics. 



1 INTRODUCTION 

Binary stars play a fundamental role in the evolution of globular 
clusters for at least two important reasons. First, the evolution of 
stars in binaries, whether in a cluster or in the galactic field, can be 
very different from the evolution of the same stars in isolation. In 
a dense environment like a globular cluster, this difference is exac- 
erbated by dynamical encounters, which affect binaries much more 
than single stars. Second, binary stars crucially affect the dynami- 
cal evolution of globular clusters, providing (through inelastic colli- 
sions) the source of energy that su p ports them again st gravothermal 
collapse jGoodman & Hull 1 19891 : bao et alj Il99ll: iFregeau et ail 



2003). In the "binary burning" phase, a cluster can remain in quasi- 
thermal equilibrium with nearly constant core density and velocity 
dispersion for many relaxation times, in a similar way to that in 
which a star can maintain itself in thermal equilibrium for many 
Kelvin-Helmholtz times by burning hydrogen in its core. The bi- 
nary fraction 1 (and the initial, primordial binary fraction in partic- 
ular), is therefore one of the most important parameters that de- 
termine the evolution of globular clusters. However, most previ- 
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1 Throughout this paper, the "binary fraction" among a particular group of 
objects is defined as the number of objects that are binaries divided by the 
total number of objects. So, for example, a primordial binary fraction of 
50% implies that 2/3 of main-sequence stars are in binaries initially. How- 
ever, a binary fraction of 50% among white dwarfs later on does not imply 
that 2/3 of white dwarfs are in binaries, as some white dwarf companions 
may be main-sequence stars. 



ous dynamical studies of globular clusters — even those including 
binaries — have neglected stellar evolution, which can significantly 
impact the properties and survival of binaries and hence the reser- 
voir of energy they provide. 

At present, there are very few direct measurements of binary 
fractions in clusters. However, even early observations showed that 
binary fractions in globu lar cluster cores are smaller than in the 
solar neig hborhood fe.g.. lCote etaill 19961) . Recent Hubble Space 
Telescope (HST) observations have provided fu rther cons traints 
on the binary fractions in many gl obular clusters iBellazzini et alj 
l2002bl : lRubenstein & BailviJ 19971) . The measured binary fractions 
in dense cluster cores are found to be very small. As an example, 
the upper limit on the core binary fraction of NGC 6397 is only 5- 
7% ICool & Boltorl2002h. On the other hand, in very sparse clus- 
ters , ]^NGC2M^^^^^^^^^), but also in some other 
"core -collapsed" clusters, like NGC 6752 iRubenstein & Bailvnl 
Il997t) . the upper limit for the binary fraction can be as high as 
~ 30%. 

For the initial binary fraction in globular clusters, there are 
of course no direct measurements. However, there are no obser- 
vational or theoretical arguments suggesting that the formation of 
binaries and hierarchical multiples in dense stellar systems should 
be significantly different from other environments like open clus- 
ters, the Galactic field, or star-forming regions. Binary frequencies 
> 50% are found in the solar neighborhoo d and in open clusters 
lHalbwachs. Mayor. Udrv , & Arenovj2003j) . T Tauri stars also hav e 
a very high binary fractioiT TKoMe^Lehrert. & ZinneckeJl200lh . 
For the range of separations between 120 and 1800 AU, their bi- 
nary fraction is comparable to that of main sequence stars in the 
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solar neighb orhood jBrandner et alll99d) . while at shorter periods 
it is higher jMeldl2003h ." Furthermore, many stars are formed in 
systems of multiplicity 3 or higher: in the field their abundance i s 
no less than 40% for inner periods ^ 10 days iTokovinirll 1 9971) . 
All this suggests that, in dense stellar systems as well, most stars 
could be formed in binary and multiple configurations. 

Most dynamical interactions in dense cluster cores tend to de- 
stroy binaries (the possible exception is tidal capture, which may 
form binaries, but turns out to play a negligible role; see § 5.2). 
Soft binaries (with orbital speeds lower than the cluster velocity 
dispersion) can be disrupted easily by any strong encounter with an- 
other passing star or binary. Even hard binaries can be destroyed in 
resonant binary-binary encounters, which typi cally eject two sin- 
gle stars and leave only one binary remaining lMfkkola|^83^or 
produce physical st ellar collisions and mergers <Baconetalll996t 
iFregeau et all2004 . 

In addition, many binary stellar evolution processes lead to 
disruptions (e.g., following a supernova explosion of one of the 
stars) or mergers (e.g., following a common envelope phase). These 
evolutionary destruction processes can also be enhanced by dy- 
namics. For example, more co mmon envelope systems form as a 
result of exchange interactions iRasio. Pfahl. & Rappaportll2000h . 
and the orbital shrinkage and the development of high eccentrici- 
ties through harden ing encounters may lead to the co alescence of 
binary components <Hillsll984LlHurlev & Sharj2003l) . 

It is therefore natural to ask whether the small binary fractions 
measured in old globular clusters today result from these many de- 
struction processes, and what the initial binary fraction must have 
been to explain the current numbers. We address these questions 
in this paper by performing calculations that combine binary star 
evolution with a treatment of dynamical interactions in dense clus- 
ter cores. In § 3 we describe in detail the method we use, follow- 
ing a brief overview of the theoretical background in § 2. We test 
our simplified dynamical model by comparing it against full Monte 
Carlo iV-body simulations in § 4.1. In § 4.2 we use semi-analytical 
estimates to predict the upper limit for the final binary fraction in 
dense clusters. In § 4.3 we estimate the lower limit for the final bi- 
nary fraction and analyse which mechanisms of binary destruction 
are most efficient as a function of cluster age. In § 5 we present our 
numerical results for the evolution of the binary fraction in dense 
cluster cores, and we compare these results with observations. In 
particular, using our theoretically predicted period distribution, we 
re-examine observations of 47 Tuc and re-derive constraints on the 
core binary fraction. In the final discussion (§ 6), we point out how 
our results may be helpful in interpreting observations of core bi- 
nary fractions in other clusters, and we discuss the required ini- 
tial conditions for simulations of clusters with binaries, as well as 
which methods are best suited for these simulations. 



2 BINARY POPULATION SYNTHESIS WITH 
DYNAMICS 

There are several possible ways to approach the study of binary 
evolution in dense clusters. The traditional approach is to start from 
full TV-body simulations to study the dynamics of the stellar system 
and introduce on top of this various simplified treatments of single 
and binary star evol ution. This has been used for many years (for 
recen t examples see IPort egies Zwartetal.ll200ll :lsh ara & Hurlevl 
12001 iHurlev & Sharall2003h . Unfortunately, even with the fastest 
special-purpose computers available today, this direct TV-body ap- 
proach remains extremely expensive computationally, so that pre- 



vious studies have been limited to small systems like open clus- 
ters and with limited coverage of parameter space. In addition, be- 
cause binaries are particularly expensive to handle computationally 
(as they increase enormously the dynamic range of direct TV-body 
simulations), these previous studies have also been performed with 
unrealistically small numbers of binaries. For example, the time re- 
quired to perform just one direct iV-body simulation of a cluster 
containing 2 x 10 stars with all stars formed initially in binaries 
would be at least a year on the GRAPE-6, with some dependence 
on the initial binary parameters 2 

Alternatively, a binary population synthesis code (e.g., 
iHurlev et all2002l) . normally used to evolve large numbers of stars 
and binaries without dynamical interactions, can be extended by 
introducing a simple treatment of dynamics. In this type of ap- 
proach it is often assumed that all the relevant parameters of 
the cluster (e.g., central density and velocity dispersion) remain 
constant throughout each dynamical simulation, i.e., the dynam- 
ics is assumed to take place in a. fixed background cluster. Many 
previous studies of dense stellar sys tems have been based on 
this type of approximation (see, e.g., iHut. Mc Millan. & Romanil 
1992tlDi Stefano & Rappaportlll994LISigurdsson & Phinnevl 1995 1 ; 
Portegies Zwart Hut. McMillan & VerbunJl997albHDaviesl 19951 
Davies & Benzl Il995t [Davies! Il997t IRasio. Pfahl. & Rappaportl 
l200Ct ISmith & Bonnellll200l[). This approach , sometimes called 
"encounter rate technique" lBenacauistal2002T) . is computationally 
much less expensive than direct iV-body simulations and hence al- 
lows the systematic exploration of the vast parameter space of ini- 
tial conditions for clusters and their primordial binary populations. 
In addition, the use of sufficiently large numbers of stars and bi- 
naries makes the simulations more realistic. Although obviously 
much less accurate in its description of the overall cluster dynam- 
ics, this method opens the possibility of studying "star cluster ecol- 
ogy" in considerably greater detail than has been possible with N- 
body simulations. In particular, it makes it possible to study in de- 
tail the rare but important evolutionary channels that may play a 
crucial role in the formation of some of the most interesting trac- 
ers of dynamical interactions in dense clusters, such as ultracom- 
pact X-ray binaries, millise cond pulsars, and cataclysmic variables 
jlvanova & Rasiol2004jbl) . 

Unfortunately, it is difficult to compare these two approaches, 
as each is based on a very different set of simplifying assumptions. 
There are no comprehensive studies of dense stellar systems includ- 
ing a self-consistent treatment of both dynamics and binary star 
evolution. In many recent TV-body simulations for large clusters 
(using either Aarseth-type codes or Henon's Monte Carlo method; 
lAarsethlEobll: IFregeau et ai]|2003l) . binary stars are treated in the 
point mass limit and soft binaries are eliminated from the start. Bi- 
nary destruction can then occur only through resonant 4-body in- 
teractions. However, iV-body studies of open clusters that incorpo- 
rate realistic treatments of binary stellar evolution have shown that 
stellar evolution affects the binaries significantly, and that, even in 
these low-density environments, the complex interplay between bi- 
nary evolution and dynamics, even for soft binaries, can play an 
important role in the overall cluste r evolution and in dete rmining 
the properties of surviving binaries iHurlev & Sharj2003l) . 

The second approach, "binary population synthesis with dy- 
namics", which we have adopted in this work, suffers from the lack 



2 J. Hurley, personal communication. The estimate is based on 5 days re- 
quired to si mulate an open cluster of 20000 stars with 2000 binaries on the 
GRAPE-6 in lShara & Hurlevl |200A 
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of self-consistent dynamical evolution of the cluster, which is as- 
sumed to remain in a constant state of thermal equilibrium for its 
entire evolution. This state, where the energy production through 
"binary burning" in the core is balanced by the outer energy flux 
into the cluster halo, does indeed provide nearly constant condi- 
tions throughout a typical globular cluster's lifetime. The excep- 
tion might be "core-collapsed" clusters, which may have run out of 
binaries and evolved to a much more centrally concentrated state. 
Typically, the density and "temperature" profiles of a cluster do 
not change much as long as there are enough binaries remaining to 
provide support against gravothermal contraction. Stellar interiors 
provide a useful analogy: as long as a star keeps enough hydrogen 
to burn in its core, it can remain in thermal equilibrium on the main 
sequence and avoid core contraction and envelope expansion. Just 
like main-sequence stars, globular clusters can maintain a nearly 
constant interior structure for many thermal time-scales (i.e., many 
relaxation times) as long as they do not run out of "fuel" (binaries). 
This behavior is expected qualitatively (see, e.g.. lGoodman & Hutl 
Il989l) . and has now been demonstrated quantitatively in many dif- 
ferent studies using various numerical t echniques for cluster simu- 
lations. For example, the recent study bv lFregeau et alj f2003l) con- 
sidered the evolution of idealized clusters of equal-mass stars with- 
out stellar evolution for a range of initial binary fractions. For the 
case of an isolated Plummer model with 10% initial hard binaries 
they found (see their Fig. 4) that the core radius of this cluster can 
remain nearly constant (to within a factor ~ 2) for many tens of 
half-mass relaxation times (i.e., more than a Hubble time for most 
Galactic globular clusters, where the half-mass relaxation time is 
~ 10 9 yr). For a more realistic cluster model with 20% hard bina- 
ries initially and tidal truncation, they found that, after about 40 t^, 
when the cluster is about to disrupt in the Galactic tidal field, the 
core radius still has not varied by more than a factor ~ 2 over the 
entire evolution (see their Fig. 11); and over the first 10 t r h, the core 
radius changed by less than ~ 20%. The central velocity dispersion 
also d oes not vary much with time (see, e.g.. lGiersz & Spurzeml 
120031 Fig. 1). Similar results have been obtained from direct N- 
body simulations of open clusters, where the central density and 
velocity dispersion also rem ain nearly constant in m odels with sig- 
nificant numbers of binaries <Hurlev& Shani2003h . 



There are several binary population synthesis codes in use to- 
day. Only a few of them include, in addition to detailed treatment 
of mass transfer phases, stellar evolution along with tidal int erac- 
tion of binary components. The code of iHurlev et alj l2002h was 
designed to study low- and intermediate-mass st ars leading to the 
formation of white dwarf systems. The code of iBelczvnski et aD 
l2002Ll2005l in preparation) was originally developed to investigate 
more massive stars — progenitors of neutron star and black hole sys- 
tems (calibrated against full binary stellar evolution codes for mass- 
transfer phases) — and was expanded recently to treat carefully bi- 
naries with white dwarfs as well. We use the latter code (called 
StarTrack) to follow the evolution of single and binary stars in 
our simulations. As mentioned earlier, we also treat all dynamical 
encounters explicitly by direct integration, using a recently devel- 
oped numerical toolkit for small- AT gravitational dynamics that is 
partic ularly well suited to performing 3-body and 4-body integra- 
tions jFregeau et all2004h . 



3 METHODS AND ASSUMPTIONS 

3.1 Cluster Initial Conditions 

Our initial conditions are described by the following parameters: 
total number of stars iV (single or in a binary), initial mass function 
(IMF), binary fraction ft, distribution of binary parameters (period, 
P, eccentricity, e, and mass ratio, q — m^/mi < 1). We typi- 
cally adopt standard choices used in previous population synthesis 
studies, which are based on available observations for stars in the 
field and in young star clusters (Sills et al. 2003). For most of the 
calculations reported here, we use the following "standard" initial 
conditions: 

• We adopt the IMF of lKroiipl <2002l) . which can be written as 

a broken power law dN oc m~ a dm, where a = 0.3 for tti/Mq < 
0.08, a = 1.3 for 0.08 < m/M Q < 0.5, a = 2.3 for tti/Mq > 
0.5. We assume that all stars are formed in a single burst of star 
formation at t = in our simulations. 

• We consider a wide range of stellar masses from 0.05 Mq to 
lOOAfo 3 

• The binary mass ratio, q, is assumed to be distributed uni- 
formly in the range < q < 1. This i s in agreement with ob- 
servations for q > 0.2 JWoitas et all200lh . 

• The binary period, P, is taken from a uniform distribution in 
log 10 P over the range P — 0.1-10 7 d. 

• The binary eccentricity, e, follows a thermal distribution with 
probability density p(e) — 2e. 

• We reject systems where binary components would overflow 
their roche lobe at pericenter. 

The initial average stellar mass is then (m) ~ 0.5 Mq, and, 
with the flat mass ratio distribution, the average binary mass is 
(m b > ^ 0.75 Mq. 

3.2 Stellar Evolution 

We evolve all stars (single and binary) using the population syn- 
thesis code StarTrack I Belczvnski et al ] |2002h . which has re- 
cently been updated significantly iBelczvnski et all2005l in prepa- 
ration). This is the only current population synthesis code that in- 
corporates detailed treatments of all possible types of mass trans- 
fer (MT) episodes: stable MT (conservative or non-conservative), 
unstable MT (thermally or dynamically), and thermal time-scale 
MT. Also included are the effects of Eddington-limited mass ac- 
cretion and transient behavior of accretion discs (based on the disc 
instability model). StarTrack allows us to follow the evolution 
of binaries with a large range of stellar masses, metallicities, and 
star formation histories (constant rate, sudden or exponential bursts, 
etc.). StarTrack also models in detail the loss of mass and an- 
gular momentum through stellar winds (dependent on metallicity) 
and gravitational radiation, asymmetric core collapse events with 
a realistic spectrum of compact object masses, and the effects of 
magnetic braking and tidal circularization on close binaries. In our 
simulati ons, we adopted the new prescription for magnetic braking 
givenbyjlvanova & TaamH2003h . Compared to the older prescrip- 
tion JVerbunt & Zwaanll98ll) closer binaries lose their angular mo- 
mentum at a slower rate and hence survive as binaries longer. The 
evolution of single stars in StarTrack is based on the analytic 



3 The lower limit is chosen in order to provide a self-consistent mass-ratio 
distribution for binaries with primaries down to 0.15 Mq. 
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fits provided bv lHurlev et al] fcOOOl) . but includes a more realisti c 
determination of compact object masses iFrver & Kalogeral200ll) . 

We treat the evolution of stellar collision and binary merger 
products following the general "rejuvenation" prescription of 
iHurlev et alJbOOa) . It ensures that the merger product has the same 
total amount of hydrogen, helium, and carbon as the two parent 
stars together. For some stars, the assumptions made in the treat- 
ment of the merger depend on the types of stars and the type of 
merger (collision vs binary coalescence). For example, we assume 
that there is no accretion on to a neutron star during a physical 
collision, and that the other star, if it is unevolved, is destroyed 
completely (e.g., we do not consider the possible formation of a 
Thorne-Zhytkow object). We treat as a "dynamical common en- 
velope" (CE) event the outcome of a physical collision between 
a compact object and a red giant, applying a standard "alpha pre - 
scription" (where we adopt o?ceA = 1; see llben & Livioll 19931) . 
but taking into account the initial positive energy. In particular, if 
this compact object is a neutron star, a compact binary contain- 
ing a neutron star and a white dwarf is formed. We assume that, 
during the CE phase, the neutron star will accrete a significant 
amoun t of the envelope mat erial and will become a millisecond 
pulsar iBethe & Brown 1998). If the resulting mass of the neutron 
star exceeds the limit for a neutron star (taken in our simulations to 
be 2A/q), we assume that a black hole is formed. 

To evolve the cluster population of single stars and binaries 
in our code, we consider two basic time-scales. One is associ- 
ated with the evolutionary changes in the stellar population, At ev , 
and the other with the rate of encounters, Ai C oli (see § 3.4). The 
evolution timestep At cv is computed so that no more than 2% of 
all stars change their properties (mass and radius) by more than 
5%. The global timestep for the cluster evolution is taken to be 

At = min[tev,*dyn] 



3.3 Dynamical Cluster Model 

As we described in § 2, our model for the cluster dynamics is 
highly simplified. We adopt a simple two-zone, core-halo model 
for the cluster. We assume that the core number density, n c , and 
one-dimensional velocity dispersion, ctid, as well as the half- 
mass relaxation time, t T h, remain strictly constant throughout the 
evolution. While dense globular clusters of interest have <7id ~ 
lOkms -1 , the core density can vary by several orders of magni- 
tude. Here we set n c = 10 5 pc -3 for most calculations, representa- 
tive of a fairly dense cluster like 47 Tuc. In general, n c is the main 
"knob" that we can turn to increase or decrease the importance of 
dynamics. Setting n c — corresponds to a traditional population 
synthesis simulation, where all binaries and single stars evolve in 
isolation after a single initial burst of star formation. To model a 
specific cluster, we match its core mass density today, p^, central 
velocity dispersion, and half-mass relaxation time. 

The escape speed from the cluster core can be estimated from 
observations as v c ~ 2.5 030 JWebbinkll 19851) . where 03D is the 
three-dimensional central velocity dispersion. Following an inter- 
action or a supernova explosion, any object that has acquired a re- 
coil speed exceeding v e is removed from the simulation. Acquiring 
a large recoil velocity in a dynamical encounter is a very efficient 
mechanism for ejecting low-mass objects from the cluster. We find 
generally that recoil to the halo does not play a significant role: the 
recoil velocity into the halo differs by ~ 10% from the escape ve- 
locity from the cluster, and affects only a small number of objects. 

For computing interactions in the core, the velocities of all 



objects are assumed to be distributed according to a lowered 
Maxwellian <KinsJl965l) . with 

f(y) —v 2 /a 2 {m) [exp( — 1.5« 2 /cr 2 (m)) — exp( — 1.5v 2 /<r 2 (m))] , (1) 

where <r(m) = ((m}c/m) 1 ^ 2 a3D (assuming energy equipartition 
in the core) and v e is the escape speed. Here (m) c is the average 
mass of an object in the core. In addition, we use a to impose a cut- 
off for soft binaries entering the core: Any binary with maximum 
orbita l spee d < 0.1<T3d is immediately broken into two single stars 
lHillsll99Ch . 

In the presence of a broad mass spectrum, the cluster core is 
always dominated by the most massive objects in the cluster, which 
tend to concentrate there via mass segregation. As stars evolve, the 
composition of the core will therefore change significantly over 
time. M ass s egregation in gl obular clusters was investigated re- 
cently in lFregeau et alj 120021) by considering both light and heavy 
tracers in two-component models. It was found that the character- 
istic mass-segregation time-scale is given by 

t sc ~ 10C((m>h/m)<rh ■ (2) 

Here C is a constant of order unity and (m)h is the average mass 
of an object in the halo, and m is the current mass of the ob- 
ject. This equation represents a diffusion process and can be ap- 
plied to all stars, not just to those more massive than average. 
Indeed, even low-mass objects may (rarely) diffuse into the clus- 
ter core on a long time-scale, although on average they will tend 
to drift outwards. However, for very light objects with masses 
< 0.4(m)h (which is typically ~ 0.3 Mq at the beginning and 
~ 0.15 Mq after ~ 10 Gyr), this expression becomes less accurate, 
although it remains valid qualitatively. For a more recent discussion 
of mass segregation in the prese nce of a broad mass spectrum, see 
iGurkan. Freitag. & Rasiolj2004h . 

To model mass segregation in our simulations, we adopted the 
time-scale given by equation 0, but treated the process as stochas- 
tic: the probability for an object of mass m to enter the core after a 
time t s is sampled from a Poisson distribution, 

p(t s ) = (l/t ac ) exp(-t a /t ac ). (3) 

This treatment ensures that all stars heavier than ~ 0A{m)h dif- 
fusing into the core will have the appropriate mass spectrum and 
that interactions will occur between objects drawn from the correct 
distribution. 

Eq. was derived for the restricted case of a two-component 
cluster - without a realistic IMF - therefore C is unknown by a 
factor of a few. We find in simulations that the final core mass is 
nearly proportional to 1/C(ra)h. In order to obtain a better fit to 
observations for the core mass versus total cluster mass relation, we 
have fine-tuned eq. using data for 47 Tuc, in particular the ratio 
of the core mass to the total mass of the cluster. For this cluster we 
adopted a core mass of 10 M r and a total cluster mass of 1 6 Mn 
at pre sent iFreire et alj|200ll ) ; we also take t r h = 10 9 yr iHarrisl 
Il99d) and an age of 11 Gyr iGratton et ai]|2003l) . While the core 
mass can be found directly from our simulations, the total cluster 
mass has an uncertainty due to the IMF cut-off at the low mass 
end in our standard cluster model. First, we found the total mass 
of a cluster model evolved to 1 1 Gyr and the initial number of very 
massive stars (defined as those producing a black hole at the end 
of their stellar evolution). We find that at 11 Gyr, the cluster has 
145Mq per black hole (or per heavy primordial star), when the 
IMF extends down to 0.01 Mq. This allows us to normalize our 
model to the real cluster mass: with this ratio we have an estimate 
for the total cluster mass when we use a higher cut-off for the IMF 
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(0.05 Mq). This now gives us the ratio of the core mass to the total 
cluster mass corresponding to our simulations. We find that the best 
fit for 47 Tuc gives C(m) h = 3M . 

3.4 Treatment of Dynamical Interactions 

All objects in our simulations are allowed to have dynamical in- 
teractions only after they have entered the cluster core. We use a 
simple Monte Carlo prescription to decide which pair of objects 
actually have an interaction during each timestep. 

The cross section for an encounter between two objects, of 
masses mi and mj, with relative velocity at infinity Vij, is com- 
puted as 

Si.j = ^d 2 max (l + Vp/vu) , (4) 

where d max is the maximum distance of closest approach that de- 
fines a significant encounter and Vp — IGijrii + mj)/dmax is the 
velocity at pericenter. Here the index i (lowercase) reflects an in- 
dividual object in the core, while J (uppercase) denotes a random 
representative object from the subclass of objects J. In order to 
more accurately determine encounter rates, at each timestep bina- 
ries and single stars in the core are divided into 100 subclasses: 10 
by size (radius for single stars or semimajor axis for binaries) and 
10 by mass. Boundaries between mass-subclasses are fixed approx- 
imately as 0.1 x 2 n . Subclasses by size depend on the current sizes 
of single and binary populations (separately), with the step between 
subclasses taken as 5 log 10 r hin = 0.1 (log 10 i? max - log 10 7? min ). 
The encounter rate for a given object i and an object from subclass 
J is T^j = njSijVij, where nj is the number density of objects 
in subclass J, and the cross section and relative velocity are defined 
for an average object in subclass J. 

The total interaction rate for a given object i is the sum of the 
interaction rates with all relevant subclasses, T, = ^ njSijVij. 
The corresponding interaction time is r, = The actual time 

for an encounter ti follows a Poisson distribution with mean n. In 
practice, we generate a random number < X < 1, and assume 
that the encounter happened if ti = Xr, ^ At. The timestep is 
limited so that Atd yn ^ 0.25 mini Ti. We keep track separately of 
the time-scales n t j for interactions with each subclass J, and the 
corresponding ti t j = X.jTi t j is generated from an independent 
random number. If an encounter happened, it is assumed to be with 
an object from the subclass with the smallest ti.j. The actual inter- 
acting object J from that subclass J is randomly selected from the 
list of non-interacted objects in this subclass. 

In this paper we consider separately binary-binary, binary- 
single, and single-single interactions. The cross sections and rates 
are calculated using d m ax = 5(Ri + (R).j) for single-single, 
d m ax = 3(6, + (R)j) for binary-single and d max = 3(6, + (&),/) 
for binary-binary. Here bi = a,(l + e,) is the apocenter separation 
of the binary (a^ is the semi-major axis and e, is the eccentricity), 
(R)j is the average stellar radius in the subclass J, and (6) j is the 
average apocenter separation in the subclass J. A single-single in- 
teraction with pericenter distance di ^ 2(Ri + R.j) is treated as a 
physical collision and assumed to lead to a merger or a dynamical 
CE phase. If 2 < di/(Ri + R.j) < 5, we check whether a binary 
co uld be produced via tidal c a pture using the approach described 
in IPortegies Zwart & Meinerl 1 1993l> . If a tidal-capture binary is 
formed, its eccentricity is set to zero and semi-major axis set to 2di, 
assuming rapid circular ization (~ 10 yr) as predicted by the sta n- 
dard model described in McMillan. McDermott. & Taaml ll987l) . 

Each dynamical interaction involving a binary is calculated 
using Fewbody, a new numerical toolkit for simulating small-TV 



gravitational dynamics that provides automatic calculation termi- 
nation and classifica tion of outcomes (for a detailed description see 
iFregeau et all2004l) . Fewbody numerically integrates the orbits of 
small-TV systems, and performs collisions in the sticky-star approx- 
imation. Fewbody's ability to automatically classify and terminate 
calculations as soon as the outcome is unambiguous makes it well- 
suited for carrying out large sets of binary interactions, for which 
calculations must be as computationally efficient as possible. 



4 TEST CALCULATIONS AND SIMPLE ESTIMATES 

4.1 Comparison with TV-body Simulations 

We have compared our simple dynamical model to recent re- 
sults from fully self-consistent Monte Carlo simulations based 
on Henon's algorithm for solving the Fokker-Planck equa- 
tion jjoshi. Rasio. & Portegies ZwarJl2000t lloshi. Nave. & Rasiol 
1200 J) . For our test we used the results obtained for an idealized 
model cluster containing 20% primordial hard binaries (binding en- 
ergies in the range 1 — 133 kT, where kT is the average kinetic 
energy of an object in the cluster) for a King model wi th dimen- 
sionless central potential Wo = 7 iFregeau et al]l2003l hereafter 
F03). In this simulation all stars had equal masses, were treated as 
point masses (no physical collisions), no stellar evolution was taken 
into account, and all binary interactions were treated using simple 
recipes. 

We used our code to perform a similar simulation: we consid- 
ered a cluster consisting of equal-mass stars, and we turned off all 
stellar evolution and physical collisions. To fit the F03 model, we 
took the core mass as 8.3% of the total cluster mass (corresponding 
to a King model with Wo = 7), and we took an average star mass 
of 1 Mq and cr r = 10 km/s. For a total cluster mass of 3 x 10 5 Mq 
these conditions imply — 8 x 10 8 yr, and n c = 2000 pc -3 . We 
evolved the cluster for 20 t r h. 

In Fig. \\\ we show our results for the core and halo binary 
fractions as a function of time, compared with the model of F03. 
The agreement with F03 is excellent: the core binary fraction rises 
very quickly to ~ 35% and then remains close to this value for 
~ 10t r h. In contrast, the halo binary fraction decreases more grad- 
ually from 20% to 10%. Considering the differences between the 
two treatments, especially for binary interactions, this agreement is 
quite remarkable. 

In the bottom panel of Fig.Qwe show the evolution of the core 
radius r c . The core radius evolves in the same way as in the dynam- 
ical simulation of F03. Overall, the core radius does not change 
much during the entire evolution and its value is consistent with 
the measured values for observed globular clusters with similar pa- 
rameters (e.g., NGC 3201 or NGC 6254, in which the total cluster 
mass, the central number density, and the central velocity disper- 
sion are similar to those in our model). 

4.2 Semi-analytic Estimates 

One can estimate the final binary fraction /b, c in a dense environ- 
ment by considering several mechanisms of binary destruction: 

• all soft binaries are usually destroyed during a strong en- 
counter. 

• some fraction of hard binaries is destroyed through stellar evo- 
lution (mergers or disruptions after supernova explosions) 

• when a hard binary has a strong encounter with another hard 
binary or a single star, it can exchange its less massive component 
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Figure 1. Evolution of the core and halo binary fractions (top) and the core 
radius (bottom) in our test model, compared with the F03 model (see text). 
In the top panel, the solid line shows the binary fraction in the core and the 
dashed line shows the binary fraction in the halo, both for the test model. 
The dotted line shows the binary fraction in the core in the F03 model and 
the dash-dotted line shows the binary fraction in the halo in F03 model. In 
the bottom panel the solid line shows the core radius in the test model and 
the dotted line shows the core radius in the F03 model. 
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Figure 2. Destruction of primordial binaries by stellar evolution, shown in 
the parameter space of total initial binary mass and initial binary period. 
Solid lines are lines of constant binary hardness and dashed lines are lines 
of constant collision time. 



for a more massive star, shrink its orbit, or be destroyed in a colli- 
sional merger. 

It should be noted that in our simplified semi-analytical treatment 
we neglect the effect of mass segregation, which tends to increase 
the core binary fraction. (This issue is discussed in more detail in 
Section^) 

Let us consider a dense environment with number density 
10 5 pc~ 3 , oid = lOkm/s and with an average mass of 0.5 Mq. 
With our choices of initial parameters, and an average mass of 
0.5Mq, 40% of all primordial binaries initially are soft (this frac- 
tion would be 50% with respect to the average mass of IMq). We 
introduce r\ - the hardness of a binary system — as 



V = 



(5) 



aa 2 (m) ' 

where a is the binary separation, m 1 and m2 are the masses of the 
binary components, and (m) is the average mass of a single star. 
Binaries that have r\ < 1 are termed soft, and those with r\ > 1 are 
termed hard. 

To find how many hard binaries will be destroyed by stellar 
evolution alone, we calculated the probability of binary destruction 
as a function of its initial total mass and orbital period, using the 
binary population synthesis code (see Fig|2j. This simulation was 
done with 1.25 x 10 5 binaries distributed initially flat in log Pd and 
log Aftot (in order to have better resolution for destruction rates 
for high mass binaries), where Pd is the binary period in days and 
Mtot = Mi + M2 is the total binary mass in Mq , and was evolved 
for 14 Gyr. 

The result is striking: most of the very hard binaries, with 
hardness ratios 77 > 100, are destroyed by stellar evolution. The 
empty space near the bottom right corner of Fig |2| reflects the ab- 
sence of systems below the minimum period for binaries on the 
main sequence (the period at which stars come into contact). For 
binaries with total mass ^ 10 M@ destructions mainly occur dur- 
ing the first 10 s years of cluster evolution. Binaries with period 
^ 10 4 days are mainly destroyed through SN explosions. Binaries 
with period ^ 10 days are destroyed mainly via mergers at the MS 
stage. For periods 10-10 4 days the destructions are associated with 
common envelope evolution and occur at later times. Destructions 
in binaries of smaller masses are not much different from more 
massive binaries at small periods, but are not destroyed through SN 
explosions at large periods. Also, the CE event in less massive bi- 
naries occurs when the donor is a less evolved giant (at the first red 
giant branch). 

One may expect that destruction rates should vary smoothly; 
however, binary evolution involves many qualitatively different 
events. In particular, an interesting fluctuation in destruction rates 
can be seen at logP ~ 1.75 and logMtot ~ 0.95, where local 
destruction rates are lower than for nearby binaries. For these and 
nearby binaries the destruction rates are about the same for mass 
ratio q < 0.5 but different for larger q. For binaries of masses close 
to and smaller than these, most destructions for q > 0.5 occur dur- 
ing the CE event between a WD and a giant. This CE phase is the 
second interaction in the binary and follows the stable MT event 
with the other donor. When the total binary mass is smaller, the 
WD mass is < 0.9 Mq, and CE event leads to a merger. For bi- 
naries of higher total mass, the second interaction occurs between 
a He star or a WD and a star at the Hertzsprung gap. This MT is 
unstable and leads to a merger. This, therefore, provides for a small 
local decrease in destruction rates. 

With our IMF and considering binaries of all masses, the 
fraction of hard binaries destroyed during their evolution is 20%. 
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Among binaries where at least one star is more massive than 
0.5 M© the destruction fraction is 50%, and for those with total 
initial mass above 1 Mq, this fraction is closer to 60%. The frac- 
tion of hard binaries that is not destroyed but is instead softened 
by evolution is very small, ^ 1%. In §4.3 we will discuss in more 
detail how binary destruction rates change with time. 

Therefore, even before any dynamical processes are consid- 
ered for the hard binary population, we estimate that the final binary 
fraction cannot be higher than about 50%, and for binaries with at 
least one star more massive than 0.5 Mq this upper limit becomes 
30%. For relatively massive binaries, with total initial mass above 
1 Mq, the upper limit for /b. c is only 24%. Overall, this estimate 
already shows that (i) the expected final binary fraction in a dense 
star cluster will be low and (ii) stellar evolution cannot be neglected 
when estimating binary fractions from dynamical models of dense 
star clusters. 

Let us now consider the effects of dynamical interactions. The 
time-scale for a binary to undergo a strong encounter with another 
single star, the collision time, can be estimated as r co ii = 1 /naVoo ■ 
Assuming that the strong encounter occurs when the distance of 
closest approach d ma x ^ ka with k ~ 2, we obtain 



Toon = 3.4 x 10 13 yr k~ 2 p- 4/3 M^VW * 



(6) 



1 + 913 



(Mtot + (Af)) 



d -"-Hot U W ' 

Here (M) is the mass of an average single star in Mq, Vio = 
Voo/(10 km/s), where is the relative velocity at infinity and 
ris = n/ (10 5 pc~ 3 ), where n is the number density. Using eq. 
eq. can be rewritten in a more convenient form: 

(M) 2 



Tcoii = 1.7 x 10 s yr rfk 2 n 5 1 



Ml Ml 



(7) 



2M^+JM> 
k M^M 2 

It can be seen from eq. that the collision time for a binary 
with <r\ = 1 and Mtot = 1 Mq is only ~ 10 s yr. Overall, with our 
IMF, ~ 50% of all hard binaries have collision times shorter than 
~ 10 Gyr (see also Fig.0, and 15% have collision times as short as 
~ 1 Gyr. In addition, most of the hard binaries that could have ex- 
perienced an encounter are binaries with 77 = 1 — 100, and therefore 
they are from the population that is destroyed by stellar evolution to 
a lesser degree than binaries harder than r\ = 100. While this esti- 
mate considered only binary-single encounters, binary-binary en- 
counters will further enhance the destruction rate. Moreover, when 
the binary fraction is higher than ~ 25%, binary-binary encounters 
dominate over binary-single encounters. 

Each binary-single encounter with a hard binary can result 
in a hardening of this binary, an exchange of a companion with a 
more massive single star, or binary destruction in a physical colli- 
sion. The probability of a physical collision d uring an encounter in - 
creases strongly as the binary becomes harder iFreseau et all 2004). 
In addition, a very hard binary can be ejected from the core if its 
recoil velocity exceeds the escape speed from the cluster. Each of 
these three processes, directly or indirectly, leads to the depletion 
of binaries (immediate or delayed): acquiring a more massive com- 
panion, as well as orbital shrinkage or the eccentricity increase, 
makes systems more likely be destroyed through stellar evolution. 
As a conservative lower limit, we assume that half of the interacting 
hard binaries will be destroyed as a result of an encounter (imme- 
diately or later). This is clearly a lower limit, as scattering experi- 
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Figure 3. Binary destruction rates as a function of time for a field popu- 
lation (i.e., without dynamical interactions). Rates are given as numbers of 
events per binary per Gyr for mergers and disruptions (following a super- 
nova explosion). Note the peak in evolutionary disruptions at t ~ 10 7 — 10 8 
yr, which is due to supernovae. 



ments show that, for hard binaries containing main- sequence stars 
most encounters will lead to a physical collision iFregeau et alj 

Hooi. 

Taking everything into account, we conclude that 64% of all 
binaries will be destroyed - only k a — 36% of primordial bina- 
ries can survive to the present. The expected final binary fraction 
is therefore f h = N h /{N S + N b ) = k s /(2 - fc„) = 22% (for 
stars in the entire mass range), and it is 13% (k s = 23%) for bi- 
naries with at least one star more massive than 0.5 Mq, while for 
Mtot > 1M0 it is only 10% (k s = 18%). We note again that 
this upper limit for the final binary fraction does not take into ac- 
count several other mechanisms leading to binary destruction, such 
as binary-binary encounters, which are more likely to cause binary 
destruction than binary-single encounters. 



4.3 Encounters in Dense Environments 

In order to verify our encounter rates, we considered the evolution 
of a binary population that is completely immersed into a dense en- 
vironment with n c = 10 5 pc~ 3 for its entire evolution. We adopted 
a velocity dispersion crirj = 10 km/s and evolved the system for 
14 Gyr. This is an even simpler model of a star cluster, where all 
objects are effectively inside the "core" at all times and the effects 
of mass segregation and diffusion between a core and a halo are 
completely ignored. It allows us to study more clearly the various 
types of dynamical interactions as separate processes. We ran a se- 
quence of simulations, with (i) only binary-single encounters, with 
no physical collisions allowed during an encounter; (ii) binary- 
single encounters with physical collisions; (iii) both binary-single 
and binary-binary encounters but still without physical collisions; 
(iv) binary-single and binary-binary encounters with physical col- 
lisions; and (v) all types of encounters allowed, including single- 
single collisions. In this last case we also eliminated from the sys- 
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Figure 4. Same as Fig.0 but now for the same binary population evolved 
in a high-density environment with 10 5 binaries per pc 3 . Here disruptions 
through binary-binary (BB) and binary-single (BS) interactions dominate 
for the first ~ 5 Gyr. 



tem any object that acquired a recoil speed exceeding the escape 
velocity v e = 2.5 o"3d. 

In § 14.21 we estimated that, with only binary-single encoun- 
ters taken into account, the final binary fraction should not exceed 
22%. This assumed that all objects in the core had a mass equal 
to the average mass 0.5 Mq. With the complete IMF, a single star 
participating in the encounter can have mass higher than 0.5 Mq. 
Also, a fraction of initially hard binaries can become soft during 
the evolution, e.g., because of mass loss. As a result, we obtain 
even lower remaining binary fractions: /b. c = 16.5% for case (i) 
and /b, c = 16.4% for case (ii). These numbers are in agreement 
with the upper limits we estimated in § 14.21 

When binary-binary encounters are taken into account, the re- 
sult is even more dramatic: we obtain /b, c = 8.0% and /t>, c = 
7.9% for cases (iii) and (iv), respectively. As expected, we see that 
binary-binary encounters further enhance the binary destruction 
rate. 

When all dynamical processes are taken into account, in 
case (v), we obtain a final binary fraction /b, c = 6.9%. While our 
semi-analytical estimate of the binary fraction provided an upper 
limit, the result obtained in this section is clearly a lower limit: in a 
real cluster, not all binaries will be exposed to the high interaction 
rates of the core at all times. However, since more massive objects 
are more likely to diffuse into the core, and since their destruction 
rate can be higher than for objects of average mass, this lower limit 
is only approximate, but we may expect it to be much closer to the 
real value than our previous, conservative upper limit. 

Let us now consider in detail which mechanisms of binary de- 
struction are most important. In Fig. [5] we show the rates of binary 
destruction (events per binary per Gyr) for a field population (with 
no interactions). Except for the interval between ~ 10 7 and ~ 10 8 
years, when black holes and neutron stars are formed, the binary 
destructions are mainly coming from mergers. The power-law be- 
havior observed for the overall destruction rate at late times is on 



the one hand imposed by the rate of orbital decay driven by mag- 
netic braking, and on the other hand depends on the evolutionary 
time-scale for the companion to become a red giant. 

In Fig. [4] we show the binary destruction rates again, but for 
the dense environment. During the first few Gyr the main destruc- 
tion mechanism is the dynamical disruption of soft binaries. The 
rate of disruption through stellar evolution is smaller than in the 
field population: some massive binaries are destroyed by dynamical 
encounters before they evolve to a supernova explosion. However, 
the rate of evolutionary mergers is about the same as for the field 
for the first ~ 10 s years. This is consistent with our estimate for 
the collision time of hard binaries: at this time, hard binaries have 
not yet been hardened significantly by dynamical encounters. After 
~ 10 s yr the rate of binary evolutionary mergers is increased com- 
pared to the field population, as the dynamical hardening of hard 
binaries has now started. The rates for binary destruction through 
mergers and physical collisions are very similar: in most cases, if a 
binary is tight enough for an evolutionary merger, the most likely 
outcome of a dynamical interaction is a physical collision. 

At about 5 Gyr, the rate of binary destruction through physi- 
cal collisions starts dominating over dynamical disruptions: at this 
stage, almost no soft binaries are left, and the hard binaries are more 
likely to lead to physical collisions during an encounter. Consider 
a binary with 77 ~ 1 at this stage. From equation Q we can find 
that a collision time of 5 Gyr corresponds typically to binary com- 
ponents of 0.1 Mq and 0.05 Mq (assuming an average mass ratio 
of 0.5 and average field mass as 0.5 Mq). This binary is clearly 
at the lowest-mass end of our IMF, and therefore by this time al- 
most all soft binaries in the simulation have been destroyed. On the 
other hand, consider some more massive binary, which has evolved 
through about half of its main-sequence lifetime, and has compo- 
nent masses of 0.8 Mq and 0.4 Mq (and still the same collision 
time of 5 Gyr). The hardness of such a binary is 77 ~ 100, corre- 
sponding to a binary separation ~ 10 Rq. For such a tight binary, 
the probability of a physical collision in an encounter is almost 
100%. The total rate of binary destructions through physical col- 
lisions become comparable to that of dynamical disruptions, as can 
be seen from Fig.0 



5 RESULTS 

5.1 Overview of Cluster Models 

Initial conditions for all our models are given in Table 1. The 
first group of models is used to cover the parameter space of ini- 
tial conditions over fairly wide ranges. Our main reference model, 
Model 1, has core density n c = 10 5 pc -3 , half-mass relaxation 
time rj r h = 10 J yr, initial binary fraction /^o = 1, central ve- 
locity dispersion am = lOkms -1 , and central escape speed 
v c — 2.5<J3D = 43kms _1 . We then consider three models with 
different central densities (D3, D4, and D6), two with different half- 
mass relaxation time (T8 and T10) and one with an initial binary 
fraction decreased to 50% (B05). The initial total number of stars 
is N = 2.5 x 10 5 for all models, except for the "47 Tuc" model 
(see below), where we used N = 5 x 10 J ). In all these models the 
cluster core was assumed to contain about 1 % of the stars initially, 
and the metallicity was fixed at Z = 0.001 (these two parameters 
have very little influence on our results). 

We performed several simulations with parameters that at- 
tempt to match those of specific globular clusters in the Galaxy 
(the bottom part of Table 1). All these clusters are classified ob- 
servationally as "non-core-collapsed", meaning that they are well 
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Table 1. Initial Conditions for All Models 



Model 


log n c 


log P°M 


log t [h 


/b,0 


o"id 


Ve 


1 


5.0 




9.0 


1.0 


10. 


43. 


D3 


3.0 




9.0 


1.0 


10. 


43. 


D4 


4.0 




9.0 


1.0 


10. 


43. 


D6 


6.0 




9.0 


1.0 


10. 


43. 


T8 


5.0 




8.0 


1.0 


10. 


43. 


T10 


5.0 




10. 


1.0 


10. 


43. 


B05 


5.0 




9.0 


0.5 


10. 


43. 


M12 


3.8 


3.5 a 


9.02 


1.0 


4.5 


19.6 


M4 


4.4 


4.1 


8.82 


1.0 


4.2 


20.3 


47Tuc 


5.3 


5.1 


9.48 


1.0 


11.5 


56.8 


NGC 6388 


5.9 


5.7 


9.08 


1.0 


18.9 


78.2 



Notations: n c is the core number density in pc — 3 (assumed fixed), p'jjj 
is the observed core mass density in Mq pc -3 , t r h is the half-mass 
relaxation time in yr, f h0 is the initial binary fraction, airj is the 1-D 
velocity dispersion in kms -1 , and v c is the escape speed in kms 

a t r h for specific cluste rs are taken frojrjHarrisll l996l) , and crio from 
IPrvor & Mevlanl <1993h . and v e from lWebbinl<Hl985l) . 



Table 2. Results for Reference Models 



Model 


log PM 


/b,c 


/o.5 


fwd 


1 


4.7 


0.095 


0.15 


0.080 


D3 


2.7 


0.265 


0.37 


0.165 


D4 


3.7 


0.170 


0.25 


0.115 


D6 


5.7 


0.035 


0.055 


0.055 


T8 


4.5 


0.11 


0.13 


0.085 


T10 


4.8 


0.055 


0.07 


0.055 


B05 


4.7 


0.072 


0.010 


0.055 



Here pm is the core mass density in Mq pc — , /j, c is the binary fraction 
in the core, /o . 5 is the binary fraction for non-degenerate stars more massive 
than 0.5 Mq, and / w( j is the binary fraction among white dwarfs. Values 
for all quantities are given at 14 Gyr. 

fitted by standard King models. These are precisely the kinds of 
clusters that, theoretically, we expect to be in the "binary burn- 
ing" thermal equilibrium state. For this set, we tried to consider 
the maximum range of dynamical parameters, while concentrating 
on clusters at relatively small distances and/or very well studied ob- 
servationally, so that current or future observations could test our 
predictions. For our most important model, 47 Tuc, we also consid- 
ered additional models in which we varied the initial binary frac- 
tion or introduced a time-dependent core density (See § 5.5). For 
all our models of specific observed globular clusters, the central 
number density was chosen in order to provide (at the actual age of 
the clust er, ~ 11 — 13 Gyr) the best fit to the observed mass den- 
sity pf} (Pryor & Mevlan 1993). Metallicities for these models are 
taken from lHarria ll99o) . 

5.2 Main Reference Model 

First we present our results for a typical dense cluster, represented 
by Model 1 . With 100% binaries initially, the final core binary frac- 
tion (at 14 Gyr), /b. c , is only 9.5%. This is strikingly low, given that 
the cluster started with all stars in binaries, and that binaries should 
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Figure 5. Evolution of the core binary fraction in Model 1 . We show sepa- 
rately the binary fractions for all objects, for non-degenerate stars with mass 
^ 0.5 Mq, for non-degenerate stars with mass $C 0.3 Mq, and for white 
dwarfs. 



concentrate more into the core through mass segregation, but it is 
expected from our estimates in § 4. Decreasing the initial binary 
fraction, /b,o, to a more reasonable (but still large) 50% reduces 
/b, c further to 7%, as shown in Model B05 (see Table|2j. The de- 
pendence of /b x on /b,o is not linear. This is mainly due to mass 
segregation: decreasing /b.o also increases the ratio of mean bi- 
nary mass to mean stellar mass in the cluster, thereby resulting in a 
higher concentration of binaries in the core. 

The majority (about 75%) of destroyed binaries were dis- 
rupted by close dynamical encounters (or, rarely, following a su- 
pernova explosion). Note that some binaries that are initially hard 
eventually become soft after undergoing significant mass loss due 
to stellar evolution. About 20% of the destroyed binaries experi- 
enced mergers, typically after significant hardening through inter- 
actions. A few percent of the binaries lost were actually not de- 
stroyed but instead were ejected from the cluster as a result of recoil 
in strong encounters. Tidal capture did not play a significant role: 
the total number of tidal-capture binaries formed during the cluster 
lifetime is less than 1% of the final number of binaries in the core. 
The total core mass during most of the cluster evolution (the last 
~ 10 t r h) did not vary by more than a factor of two. 

While the final core binary fraction is extremely low, the over- 
all cluster binary fraction, which takes into account all halo bina- 
ries, remains high. However, the halo binaries consist mainly of 
very low-mass systems: the average primary mass among binaries 
remaining outside the core at 14 Gyr is 0.2 Mq, with the aver- 
age companion mass around 0.1 Mq. These binaries would be ex- 
tremely faint and hard to detect observationally. 

We have also examined different stellar subpopulations in 
the cluster: (i) the subpopulation of non-degenerate objects with 
mass (for a single star or for the primary in a binary) ^ 0.5 Mq 
and (ii) the subpopulation of all white dwarfs, single or in bina- 
ries. The binaries in group (i) may be easier to detect, e.g., from 
broadening of the main sequence in a colour-magnitude diagram 
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Figure 6. Evolution of the core binary fraction in Model 1 compared to 
models with different core number densities (D3, D4 and D6). 

iRubenstein & Bailvnll997f) . Binaries in group (i) are harder than 
less massive average binaries, so fewer of them are destroyed. On 
the other hand, stellar evolution plays a more significant role in 
the destruction of white-dwarf binaries, which were initially more 
massive and harder (Fig.|5J- Therefore the binary destruction rate 
in group (ii) is much higher, enhanced both by stellar evolution 
(mass loss and mass transfer at more advanced evolutionary stages, 
and supernovae in binaries), and by dynamical interactions (larger 
cross-section for encounters). Note also that the overall /b, c is de- 
creased partially through a lower binary fraction for degenerate ob- 
jects. 

5.3 Different Densities and Relaxation Times 

Let us now compare results for different central densities, in Mod- 
els 1, D3, D4 and D6. The evolution of /b, c for these models is 
shown in Figure [6] As expected, the core binary fraction decreases 
as n c increases. The dependence is steeper at high densities, where 
dynamical interaction rates play a more dominant role compared to 
stellar evolution. 

With respect to the half-mass relaxation time (Models T8 and 
T10), we find that, surprisingly, the model with shorter relaxation 
time has a higher core binary fraction. There are two competing 
mechanisms that play a role here: mass segregation, which brings 
binaries into the core, and dynamical interactions, which destroy 
binaries in the core. A shorter segregation time increases the rate at 
which binaries enter the core but also allows less massive binaries 
to interact. Therefore, the average mass of a binary in the core and 
the effective mass density in the core is smaller. As a result, fewer 
binaries are destroyed. 

5.4 Binary Period Distribution 

Through dynamical interactions, we would expect that the initial 
period distribution of binaries should be depleted above the bound- 
ary between hard and soft binaries. Stellar evolution should deplete 
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Figure 7. The binary period distributions in models with different densities 
(from the top: D3, D4, Model 1 and D6). Nj, is the total number of binaries 
and Nt,,p is the number of binaries in each period bin. 



Table 3. Results for Models of Specific Clusters 



Model 


Age 


log PM 


fb,c 


fo.5 


/wd 


M12 


12 


3.5 


0.170 


0.26 


0.130 


M4 


13 


4.4 


0.115 


0.175 


0.10 


47Tuc 


11 


5.1 


0.080 


0.125 


0.075 


NGC 6388 


12 


5.7 


0.06 


0.08 


0.045 



With the same notations as before, the columns give: the name of the clus- 
ter; the cluster age in Gyr; the central mass density in Mg pc - ^ in sim- 
ulations at the given age, the core binary fraction, the binary fraction for 
non-degenerate objects more massive than 0.5 Mg, and the binary fraction 
among white dwarfs. 

a fraction of hard binaries, especially at very short periods, and dy- 
namical encounters should further deplete some of the wider hard 
binaries. Indeed, for lower-density clusters, we find that the distri- 
bution remains much flatter in log P (see Fig. 0, with more and 
more of the wider hard binaries disappearing as the density in- 
creases. 

This period distribution can be used to better extrapolate ob- 
served binary fractions, which are usually limited to rather narrow 
ranges of masses and periods. In particular, it is clear that the usual 
assumption of a flat di stribution in log P fo r hard binaries at present 
in a cluster core fe.g jAlbrow et aj200lh can be very misleading. 
This will be done in § 5.5 for our models of several real clusters. 

5.5 Comparison with Observations 

We performed several simulations with parameters that attempt to 
match those of specific globular clusters in the Galaxy (Tables Q 
and [3). For these models, the total core masses we compute at 
present match well those of King models fitted to the observed 
cluster parameters, and the corresponding core radii of our models 
(from total core mass and density) are all ~ 0.5 pc, matching the 
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All binaries in the core 
MS-MS with M,,M 2 >2.5 M„ 
WD-binaries 




MS MS: Mj.Mg iO.25 M„ 
Dynamically formed 




Figure 8. Period distributions for different binaries in Model 1 (at 14 Gyr). 
The middle plot shows the period distribution of binaries containing two 
MS stars with masses greater than 0.25 Mq, bottom plot shows the period 
distribution of binaries with at least one WD; shaded area shows the fraction 
of dynamically formed binaries. Nj, is the total number of binaries and Nv, p 
is the number of binaries in each period bin. 




All binaries in the core 
MS-MS with M P M Z >0.25 M; 
WD-binaries 



MS MS: Mj.Mj, iO.25 M Q 
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WD-binaries in the core 
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Figure 9. The binary period distribution in our 47 Tuc Model, at 1 1 Gyr. The 
middle plot shows the period distribution of binaries containing two main- 
sequence stars with masses greater than 0.25 Mq; the bottom plot shows 
the period distribution of binaries with at least one white-dwarf component 
(the shaded area shows the fraction of dynamically formed binaries). As 
before Nv, is the total number of binaries and N 
in each period bin. 



b. P 



is the number of binaries 



observed values. In all cases the initial binary fraction is assumed 
to be 100%, so our results for final core binary fractions represent 
upper limits. As in all reference models, we predict low values for 
/b, c in all globular clusters, with smaller values in higher-density 
cores. 

We c ompare our derived bi nary fraction for M4 with the re- 
sults from lCote & Fischer! Il996l) . Using a Monte Carlo modelling 
of the binary population, they found that the best fit to the observed 
binaries (in the period range from 2 days to 3 years) is /b, c — 15%. 
This is in good agreement with our overall predicted /b, c = 11.5%, 
and also with our prediction for heavier main-sequence binaries, 
/o.s = 17.5%. 

In the 47 Tuc model (for which we increased the total number 
of stars initially to a more realistic N = 5 x 10 ), /b,c is only 8% at 
an age of 1 1 Gyr and about 7% at 14 Gyr (we provide this value at 
different ages given the uncertainty in observationally determined 
values for the ages of 47 Tuc; see, e.g., Schiavon et al. 2002; Zoccali 
at al. 2001; also note that /b. c does not change much over the last 
several Gyr). At first glance, this may see m to conflict with obser- 
vations. In particular. lAlbrow et alj l200ll) derive an overall binary 
fraction for the core of 47 Tuc of about 13%, from observations 
of eclipsing binaries with periods in the range P ~ 4-16 d. This 
estimate was based on an extrapolation assuming a period distribu- 
tion flat in log P from about 2.5 d to 50 yr (soft primordial binaries 
with P > 50 yr are assumed destroyed and short-period primordial 
binaries with P ^ 4 d are assumed to evolve toward much shorter 
periods through angular momentum loss by magnetic braking). In 
Figure 2 we show the final period distribution of core binaries in 
our simulation. Note that the period range of eclipsing binaries is 
near the peak of the distribution, while for longer periods the num- 
ber of binaries drops rapidly. In particular, if we concentrate on 
the binaries consistent with the observed set, with primary masses 
Mi > 0.25M© and q > 0.3, we find that the number of systems 
with periods in the range 16 d to 50 yr is about 6 times smaller than 
would be predicted by adopting a distribution flat in log P. For the 
whole period range from 2.5 d to 50 yr we have 2.2 times fewer 
binary systems compared to the flat distribution. If we take into ac- 
count this depletion of wider binaries when modelling the number 
of observed eclipsing binaries in 47 Tuc, we are led t o revise the 
observed core binary fraction from lAlbrow et al]i200ll) to 6% ± 2, 
which is much closer to our theoretical prediction. 

We performed three additional simulations for 47 Tuc, with 
fb.o = 0.25, 0.5, and 0.75. The corresponding core binary fraction 
extrapolated from observations (corrected as above) does not vary 
much among these different models. Its maximum value is obtained 
for /b,o = 0.5, which gives a core binary fraction of about 8%, and 
in this case the total number of binary systems is 1.6 times less than 
with an assumed flat distribution. 

An alternative estimate of the binary fraction in the core of 
47 Tu c is based on observations of BY Dra stars lAlbrow et alj 
l200ll) . Their estimated core binary fraction, which can be consid- 
ered a lower limit, is approximately 0.8%, 18 times lower than the 
estimate based on eclipsing binaries. This estimate was based on 
31 BY Dra binaries (observed in the period range 0.4-10 d) and 5 
eclipsing binaries (period range 4— 16 d). We analysed the core bi- 
nary population in our model in order to identify BY Dra systems 
and eclipsing binaries, and considering the ratio of the two. We 
adopted the standard definition for a BY Dra binary; prim ary mass 
in the range 0.3-0.7 M© (see, e.g.. lBopp & Fekel|[l977i) and pe- 
riod in the range 0.4-10 d (as for the observed sample in 47 Tuc). 
For eclipsing binaries we adopted the observed period range 4— 
16 d with each star 0.25A'/q. The ratio between the number of 
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BY Dra systems and the number of eclipsing binaries is found to 
be 5.9, 6.7, 7.2 and 3.5 for models with / b ,o = 1.0, 0.75, 0.5 and 
0.25, respectively. Therefore a large initial binary fraction > 75% 
is most likely. 

Of the quantities that we explicitly assume in our cluster 
model to be constant, the central (core) density is the one with the 
largest dynamic range in models that provide for dynamical evo- 
lution. Hence it is also the quantity that is most likely to signifi- 
cantly affect our final results. To test the sensitivity of our results 
to a changing central density, we have run our model of 47 Tuc 
with a central density assumed to increase by a factor of 10 from 
t = to the present. Specifically, we still match the currently ob- 
served core density while ramping up the value either exponentially 
or linearly in time, starting with a value 10 times smaller than at 
present. This could represent qualitatively the gradual core contrac- 
tion observed in some iV-body models of star clusters with binaries 
(see lAarseth & Heggidl 19931 and, for a steeper core contraction, 
see Giersz & Spurzem 2003). These models predict a present-day 
core binary fraction that is only slightly higher, by about 1-2%. To 
increase the binary fraction more significantly, the time-averaged 
core density in the cluster would have had to be many orders of 
magnitude lower than what is observed today. There is no reason- 
able dynamical history that would produce such an unusual result. 
In contrast, recent iV-body simulations show that the presence of 
just 10% hard primordial binaries leads to core radius expansion 
and therefore the core density might be higher in the past (Wilkin- 
son et al. 2003). Though we did not run another 47 Tuc model 
with a central density decreasing with time, we can predict that the 
present-day binary fraction would then be smaller by about 1-2%. 



6 CONCLUSIONS 

We have considered in detail processes of binary destruction (and 
formation) in dense stellar systems. In particular, we have shown 
that the present binary fraction in cluster cores should be relatively 
small (< 10%). This is caused not only by dynamical encounters, 
but also binary stellar evolution, which is the dominant mechanism 
for the destruction of hard binaries. We also find that the destruction 
rate due to stellar evolution is enhanced significantly by dynamical 
hardening of binaries. 

We have shown that values of binary fractions for stars in dif- 
ferent mass ranges and at different evolutionary stages can differ 
significantly. The fraction of dynamically formed binaries is higher 
when one considers stars at more advanced evolutionary stages. 
The binary period distribution evolves from flat in log P for loose 
clusters toward a sharply peaked distribution in denser clusters, 
even if all clusters have identical velocity dispersions and there- 
fore identical hard-soft binary boundaries. This implies that a flat 
period distribution should not be assumed when deriving overall 
binary fractions by extrapolation from the distribution of observed 
binaries in a narrow period range (e.g., eclipsing binaries). 

We considered several models that attempted to match ob- 
served globular clusters. For those with good available data on bi- 
naries (M4 and 47 Tuc), we found our predicted binary fractions to 
be in good agreement with observations once we take into account 
the correct binary period distribution. The main conclusion we de- 
rive from these calculations and our semi-analytic estimates is that 
the currently observed binary fractions in cluster cores suggest very 
high initial (primordial) binary fractions — close to 100%. 

In addition to their implications for the interpretation of ob- 
served binary fractions in cluster cores, our results have important 



consequences for the theoretical modeling of globular clusters us- 
ing iV-body simulations. Indeed, it is clear that "realistic" dynam- 
ical simulations of globular cluster evolution should include large 
populations of primordial binaries, with initial binary fractions in 
the range ~ 5 0%-100% (similar to what is usually assumed for the 
field; see, e.g.. lDuquennov & Mavorll99ll) . This poses a particular 
challenge for direct iV-body simulations, where the treatment of 
even relatively small numbers of binaries can add enormous com- 
putational costs. For this reason, current direct TV-body simulations 
of star clusters with large initial binary fractions typically have N 
too small to be considered representative of globular clusters (see, 
e.g., Portegies Zwart et al. 2003; Wilkinson et al. 2003). Other 
methods, such as Monte Carlo simulations, do not suffer from the 
same limitations, and routinely simulate clusters with reasonably 
large N (~ 10 5 - 10 6 ) and binary fractions (~ 30%), but have not 
yet incl uded advanced treatments of binary star evolution (see, e.g., 
Fregeau et alj2003l and references therein). 
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APPENDIX A: 
FORMATION 



ROLE OF THREE-BODY BINARY 



Consider the rate of three-body binary formation (via the close ap- 
proach of three single stars) in a dense cluster core. We denote by 
r(i5b) the rate per star of the formation of a binary with binding 
energy Eh- First, we consider T(b ^ 6 m ax) - the rate (per star) 
at which three objects come together within a region of size 6 max . 
The probability that a third object will be in the vicinity b of two 
other interacting objects is the product of the probability of the first 
two bodies meeting and the probability that during this time a third 
object will be in the same vicinity. The rate of two-body encounters 
for masses mi and ?ri2, with number density n,2 is 



T2(6 < 6max) = 7T& n 

Here 

2 2G(mi+7Tl2) 
Vr, = : 



1 + 



< V12 >' 



(Al) 



(A2) 



(A3) 



is the velocity at closest approach and 

^2 / 2 . 2n 2 „ mi + m 2 

< V12 > — (o"i + 0"2j = a < m > , 

mim2 

where a is the three-dimensional velocity dispersion. 

We define the minimum hardness for the binary to be formed 

as 

Gm\_m,2 



< m > cH 



Then 



(A4) 



(A5) 



The second object spends in the vicinity b of the first object 
a time At ~ ^ . The probability that a third object will be within 
the same vicinity is then 



p 3 ~ n 3 (b max Aiu 3 + b max ) = n 3 b n 



2^ 



(A6) 



Here V3 is the relative velocity of the third object with respect to 
the centre of mass of the first two. Effectively, it reflects velocity 
at which the population at the vicinity b is increasing. We do not 
consider here the population decrease, assuming that if the object 
was in the neighborhood, a three-body encounter has happened. 

The final result then takes the form given in 
iBinnev & Tremainel 1 19871) . §8), but with an additional mass- and 
hardness-dependent factor /, 



r 3 (7? > TJmln) 



niG < m 



/(mi,m2,m 3 ,r?) 



(A7) 
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/(mi,m 2 ,m 3 ,T)) = 

- ,2 T <S < S l- 5 (l + 2v ) (A8) 

( 1 _1_ ia.n -0 ' 5 r) m 1 m 2 

It is clear that the rate is highly dependent on the masses of the 
participating stars and decreases steeply with increasing hardness 
of the binary that is formed. 

If all encounters resulted in the formation of a binary with 
hardness rj, then the rate of binary formation would be completely 
described by equation lA7l 

Using the equation )A7t . the time-scale for three-body binary 
formation can also be written as 

T 3b = 2.08 XIO 18 ^] (C) 2 (^) 5 (iri^) 9 

ric_ nc / <m> \ 5 /<m> \ 5 o- (A9^ 
n 2 n 3 V mi / \ m 2 / i>i2 , 

r ? 5 (l + 2r;)- 1 f 1 + Ha^-O-S £~ "i^ 

Even in the core of a large, very dense cluster, with den- 
sity ~ 10 5 pc -3 and containing ~ 10 stars, no binaries with 
rj > 1 will be formed from ~ 1 Mq stars in a Hubble time. It 
has also been shown that many three-body binary formation events 
will lead to physical coll isions for stars as small as white dwarfs 
IChernoff & Huandl 199 A . We therefore neglect all three-body in- 
teractions in our cluster simulations. 



